{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a8b892c8",
   "metadata": {},
   "source": [
    "Printed and electronic copies of *Modeling and Simulation in Python* are available from [No Starch Press](https://nostarch.com/modeling-and-simulation-python) and [Bookshop.org](https://bookshop.org/p/books/modeling-and-simulation-in-python-allen-b-downey/17836697?ean=9781718502161) and [Amazon](https://amzn.to/3y9UxNb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "excited-advance",
   "metadata": {},
   "source": [
    "# Drag"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imported-table",
   "metadata": {
    "tags": []
   },
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "formal-context",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://raw.githubusercontent.com/AllenDowney/' +\n",
    "         'ModSimPy/master/modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "progressive-typing",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "plastic-trigger",
   "metadata": {
    "tags": []
   },
   "source": [
    "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. \n",
    "Click here to access the notebooks: <https://allendowney.github.io/ModSimPy/>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "superb-february",
   "metadata": {},
   "source": [
    "In the previous chapter we simulated a penny falling in a vacuum, that\n",
    "is, without air resistance. But the computational framework we used is\n",
    "very general; it is easy to add additional forces, including drag.\n",
    "\n",
    "In this chapter, I present a model of drag force and add it to the\n",
    "simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interesting-freeware",
   "metadata": {},
   "source": [
    "## Drag Force\n",
    "\n",
    "As an object moves through a fluid, like air, the object applies force\n",
    "to the air and, in accordance with Newton's third law of motion, the air applies an equal and opposite force to the object (see\n",
    "<http://modsimpy.com/newton>).\n",
    "\n",
    "The direction of this *drag force* is opposite the direction of\n",
    "travel, and its magnitude is given by the drag equation (see\n",
    "<http://modsimpy.com/drageq>): \n",
    "\n",
    "$$F_d = \\frac{1}{2}~\\rho~v^2~C_d~A$$\n",
    "\n",
    "where\n",
    "\n",
    "-   $F_d$ is force due to drag, in newtons (N), which are the SI units of force. A newton is 1 kg m/s$^2$.\n",
    "\n",
    "-   $\\rho$ is the density of the fluid in kg/m$^3$.\n",
    "\n",
    "-   $v$ is the magnitude of velocity in m/s.\n",
    "\n",
    "-   $A$ is the *reference area* of the object, in m$^2$. In this\n",
    "    context, the reference area is the projected frontal area, that is, the visible area of the object as seen from a point on its line of\n",
    "    travel (and far away).\n",
    "\n",
    "-   $C_d$ is the *drag coefficient*, a dimensionless quantity that\n",
    "    depends on the shape of the object (including length but not frontal area), its surface properties, and how it interacts with the fluid."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unable-scheduling",
   "metadata": {},
   "source": [
    "For objects moving at moderate speeds through air, typical drag\n",
    "coefficients are between 0.1 and 1.0, with blunt objects at the high end of the range and streamlined objects at the low end (see\n",
    "<http://modsimpy.com/dragco>).\n",
    "\n",
    "For simple geometric objects we can sometimes guess the drag coefficient with reasonable accuracy; for more complex objects we usually have to take measurements and estimate $C_d$ from data.\n",
    "\n",
    "Of course, the drag equation is itself a model, based on the assumption that $C_d$ does not depend on the other terms in the equation: density, velocity, and area. For objects moving in air at moderate speeds (below 45 mph or 20 m/s), this model might be good enough, but we will revisit this assumption in the next chapter.\n",
    "\n",
    "For the falling penny, we can use measurements to estimate $C_d$. In\n",
    "particular, we can measure *terminal velocity*, $v_{term}$, which is\n",
    "the speed where drag force equals force due to gravity:\n",
    "\n",
    "$$\\frac{1}{2}~\\rho~v_{term}^2~C_d~A = m g$$ \n",
    "\n",
    "where $m$ is the mass of the object and $g$ is acceleration due to gravity. Solving this equation for\n",
    "$C_d$ yields: \n",
    "\n",
    "$$C_d = \\frac{2~m g}{\\rho~v_{term}^2~A}$$ \n",
    "\n",
    "According to *Mythbusters*, the terminal velocity of a penny is between 35 and 65 mph (see <http://modsimpy.com/mythbust>). Using the low end of their range, 40 mph or about 18 m/s, the estimated value of $C_d$ is 0.44, which is close to the drag coefficient of a smooth sphere.\n",
    "\n",
    "Now we are ready to add air resistance to the model. But first I want to introduce one more computational tool, the `Params` object."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "human-cloud",
   "metadata": {},
   "source": [
    "## The Params Object\n",
    "\n",
    "As the number of system parameters increases, and as we need to do more work to compute them, we will find it useful to define a `Params` object to contain the quantities we need to make a `System` object. `Params` objects are similar to `System` objects, and we initialize them the same way.\n",
    "\n",
    "Here's the `Params` object for the falling penny:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "established-knitting",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = Params(\n",
    "    mass = 0.0025,      # kg\n",
    "    diameter = 0.019,   # m\n",
    "    rho = 1.2,          # kg/m**3\n",
    "    g = 9.8,            # m/s**2\n",
    "    v_init = 0,         # m / s\n",
    "    v_term = 18,        # m / s\n",
    "    height = 381,       # m\n",
    "    t_end = 30,         # s\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "instrumental-gross",
   "metadata": {},
   "source": [
    "The mass and diameter are from <http://modsimpy.com/penny>. The density\n",
    "of air depends on temperature, barometric pressure (which depends on\n",
    "altitude), humidity, and composition (see <http://modsimpy.com/density>). \n",
    "I chose a value that might be typical in New York City at 20 °C.\n",
    "\n",
    "Here's a version of `make_system` that takes the `Params` object and computes the inital state, `init`, the area, and the coefficient of drag.\n",
    "Then it returns a `System` object with the quantities we'll need for the simulation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "published-jesus",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import pi\n",
    "\n",
    "def make_system(params):\n",
    "    init = State(y=params.height, v=params.v_init)\n",
    "\n",
    "    area = pi * (params.diameter/2)**2\n",
    "\n",
    "    C_d = (2 * params.mass * params.g / \n",
    "           (params.rho * area * params.v_term**2))\n",
    "\n",
    "    return System(init=init,\n",
    "                  area=area,\n",
    "                  C_d=C_d,\n",
    "                  mass=params.mass,\n",
    "                  rho=params.rho,\n",
    "                  g=params.g,\n",
    "                  t_end=params.t_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "personalized-kruger",
   "metadata": {},
   "source": [
    "And here's how we call it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "executive-protection",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "portable-carbon",
   "metadata": {},
   "source": [
    "Based on the mass and diameter of the penny, the density of air, and acceleration due to gravity, and the observed terminal velocity, we estimate that the coefficient of drag is about 0.44."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "quick-cedar",
   "metadata": {},
   "outputs": [],
   "source": [
    "system.C_d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "inner-telescope",
   "metadata": {},
   "source": [
    "It might not be obvious why it is useful to create a `Params` object just to create a `System` object.\n",
    "In fact, if we run only one simulation, it might not be useful.  But it helps when we want to change or sweep the parameters.\n",
    "\n",
    "For example, suppose we learn that the terminal velocity of a penny is actually closer to 20 m/s.\n",
    "We can make a `Params` object with the new value, and a corresponding `System` object, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "great-crime",
   "metadata": {},
   "outputs": [],
   "source": [
    "params2 = params.set(v_term=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "outstanding-truth",
   "metadata": {},
   "source": [
    "The result from `set` is a new `Params` object that is identical to the original except for the given value of `v_term`. \n",
    "If we pass `params2` to `make_system`, we see that it computes a different value of `C_d`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "independent-trace",
   "metadata": {},
   "outputs": [],
   "source": [
    "system2 = make_system(params2)\n",
    "system2.C_d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "charged-explorer",
   "metadata": {},
   "source": [
    "If the terminal velocity of the penny is 20 m/s, rather than 18 m/s, that implies that the coefficient of drag is 0.36, rather than 0.44.\n",
    "And that makes sense, since lower drag implies faster terminal velocity.\n",
    "\n",
    "Using `Params` objects to make `System` objects helps make sure that relationships like this are consistent.  And since we are always making new objects, rather than modifying existing objects, we are less likely to make a mistake."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "primary-advocate",
   "metadata": {},
   "source": [
    "## Simulating the Penny Drop\n",
    "\n",
    "Now let's get to the simulation.  Here's a version of the slope function that includes drag:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "noble-stick",
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(t, state, system):\n",
    "    y, v = state\n",
    "    rho, C_d, area = system.rho, system.C_d, system.area\n",
    "    mass, g = system.mass, system.g\n",
    "    \n",
    "    f_drag = rho * v**2 * C_d * area / 2\n",
    "    a_drag = f_drag / mass\n",
    "    \n",
    "    dydt = v\n",
    "    dvdt = -g + a_drag\n",
    "    \n",
    "    return dydt, dvdt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baking-class",
   "metadata": {},
   "source": [
    "As usual, the parameters of the slope function are a time stamp, a `State` object, and a `System` object. \n",
    "We don't use `t` in this example, but we can't leave it out because when `run_solve_ivp` calls the slope function, it always provides the same arguments, whether they are needed or not.\n",
    "\n",
    "`f_drag` is force due to drag, based on the drag equation. `a_drag` is\n",
    "acceleration due to drag, based on Newton's second law.\n",
    "\n",
    "To compute total acceleration, we add accelerations due to gravity and\n",
    "drag. \n",
    "`g` is negated because it is in the direction of decreasing `y`; `a_drag` is positive because it is in the direction of increasing\n",
    "`y`. \n",
    "In the next chapter we will use `Vector` objects to keep track of\n",
    "the direction of forces and add them up in a less error-prone way.\n",
    "\n",
    "As usual, let's test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "velvet-tunisia",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_func(0, system.init, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "psychological-style",
   "metadata": {},
   "source": [
    "Because the initial velocity is 0, so is the drag force, so the initial acceleration is still `g`. \n",
    "\n",
    "To stop the simulation when the penny hits the sidewalk, we'll use the\n",
    "event function from the previous chapter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "practical-nowhere",
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func(t, state, system):\n",
    "    y, v = state\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "executive-there",
   "metadata": {},
   "source": [
    "Now we can run the simulation like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "liberal-dictionary",
   "metadata": {},
   "outputs": [],
   "source": [
    "results, details = run_solve_ivp(system, slope_func,\n",
    "                                 events=event_func)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incident-paradise",
   "metadata": {},
   "source": [
    "Here are the last few time steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "coastal-anthropology",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "suburban-martial",
   "metadata": {},
   "source": [
    "The final height is close to 0, as expected.\n",
    "\n",
    "Interestingly, the final velocity is not exactly terminal velocity, which is a reminder that the simulation results are only approximate.\n",
    "\n",
    "We can get the flight time from `results`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "interim-underground",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_sidewalk = results.index[-1]\n",
    "t_sidewalk"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "institutional-colors",
   "metadata": {},
   "source": [
    "With air resistance, it takes about 22 seconds for the penny to reach the sidewalk.\n",
    "\n",
    "Here's a plot of position as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "small-franchise",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_position(results):\n",
    "    results.y.plot()\n",
    "        \n",
    "    decorate(xlabel='Time (s)',\n",
    "         ylabel='Position (m)')\n",
    "    \n",
    "plot_position(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "plain-phone",
   "metadata": {},
   "source": [
    "And velocity as a function of time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "analyzed-criticism",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_velocity(results):\n",
    "\n",
    "    results.v.plot(color='C1', label='v')\n",
    "        \n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Velocity (m/s)')\n",
    "    \n",
    "plot_velocity(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "careful-causing",
   "metadata": {},
   "source": [
    "From an initial velocity of 0, the penny accelerates downward until it reaches terminal velocity; after that, velocity is constant."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "inside-confidence",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This chapter presents a model of drag force, which we use to estimate the coefficient of drag for a penny, and then simulate, one more time, dropping a penny from the Empire State Building.\n",
    "\n",
    "In the next chapter we'll move from one dimension to two, simulating the flight of a baseball.\n",
    "But first you might want to work on these exercises."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "planned-endorsement",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. \n",
    "You can access the notebooks at <https://allendowney.github.io/ModSimPy/>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "respective-address",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    " Run the simulation with a downward initial velocity that exceeds the penny's terminal velocity.\n",
    "\n",
    "What do you expect to happen?  Plot velocity and position as a function of time, and see if they are consistent with your prediction.\n",
    "\n",
    "Hint: Use `params.set` to make a new `Params` object with a different initial velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "inclusive-twenty",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "vertical-judge",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "greenhouse-madagascar",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "sudden-details",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "opening-jurisdiction",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "smaller-millennium",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    " Suppose we drop a quarter from the Empire State Building and find that its flight time is 19.1 seconds.  Use this measurement to estimate terminal velocity and coefficient of drag.\n",
    "\n",
    "You can get the relevant dimensions of a quarter from <https://en.wikipedia.org/wiki/Quarter_(United_States_coin)>.\n",
    "\n",
    "1. Create a `Params` object with new values of `mass` and `diameter`. We don't know `v_term`, so we'll start with the initial guess 18 m/s.\n",
    "\n",
    "2. Use `make_system` to create a `System` object.  \n",
    "\n",
    "3. Call `run_solve_ivp` to simulate the system.  How does the flight time of the simulation compare to the measurement?\n",
    "\n",
    "4. Try a few different values of `v_term` and see if you can get the simulated flight time close to 19.1 seconds.\n",
    "\n",
    "5. Optionally, write an error function and use `root_scalar` to improve your estimate.\n",
    "\n",
    "6. Use your best estimate of `v_term` to compute `C_d`.\n",
    "\n",
    "Note: I fabricated the \"observed\" flight time, so don't take the results of this exercise too seriously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "compact-bunny",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "shared-contrary",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "portable-account",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "arabic-shareware",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "valued-literature",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "sufficient-retail",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "comparable-lounge",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "ready-people",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "miniature-remark",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "frozen-termination",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
